Physical Aspects of Axonemal Beating and Swimming 



Sebastien Camalet and Frank Jiilicher 

PhysicoChimie Curie, UMR CNRS/IC 168, 26 rue d'Ulm, 75248 Pans Cedex 05, France 



O 
O 
O 



o 

CO 



43 ■ 

Oh' 
i 

o 



C/3 

o 

>v 
43 



We discuss a two-dimensional model for the dynamics of 
axonemal deformations driven by internally generated forces 
of molecular motors. Our model consists of an elastic filament 
pair connected by active elements. We derive the dynamic 
equations for this system in presence of internal forces. In 
the limit of small deformations, a perturbative approach al- 
lows us to calculate filament shapes and the tension profile. 
We demonstrate that periodic filament motion can be gener- 
ated via a self-organization of elastic filaments and molecular 
motors. Oscillatory motion and the propagation of bending 
waves can occur for an initially non-moving state via an in- 
stability termed Hopf bifurcation. Close to this instability, 
the behavior of the system is shown to be independent of 
microscopic details of the axoneme and the force-generating 
mechanism. The oscillation frequency however does depend 
on properties of the molecular motors. We calculate the oscil- 
lation frequency at the bifurcation point and show that a large 
frequency range is accessible by varying the axonemal length 
between 1 and 50pim. We calculate the velocity of swimming 
of a flagellum and discuss the effects of boundary conditions 
and externally applied forces on the axonemal oscillations. 



In addition, it contains a large number of other proteins 
such as nexin which provide elastic links between micro- 
tubule doublets, see Fig. [I]. The axoneme is inherently 
active. A large number of dynein molecular motors are 
located in two rows between neighboring microtubules 
and can induce forces and local displacements between 
adjacent microtubules (Albert et al. 1994). 

Axonemal fiagella can generate periodic waving or 
beating patterns of motion. In the case of sperm for ex- 
ample, a bending wave of the flagellum propagates from 
the head which contains the chromosomes towards the 
tail. In a viscous environment, the surrounding fluid is 
set in motion and hydrodynamic forces act on the fila- 
ment. For typical values of frequencies and length scales 
given by the size of a flagellum, the Reynolds number 
is small and inertia terms in the fluid hydrodynamics 
can be neglected (Taylor 1951). Therefore, only friction 
forces resulting from solvent viscosity can contribute to 
propulsion. Under such conditions self-propulsion is pos- 
sible if a wave propagates towards one end, a situation 
which breaks time-reversal invariance of the sequence of 
deformations of the flagellum (Purcell 1977). 
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I. INTRODUCTION 

Many small organisms and cells swim in a viscous en- 
vironment using the active motion of cilia and fiagella. 
These are hair-like appendages of the cell which can un- 
dergo periodic motion and use hydrodynamic friction to 
induce cellular self-propulsion (Bray 1992). In this paper, 
we are interested in those fiagella and cilia which con- 
tain force generating elements integrated along the whole 
length of the elastic filamentous structure. They repre- 
sent rod-like elastic structures which move and bend as a 
result of internal stresses. Examples for these systems are 
Paramecium which has a large number of cilia on its sur- 
face; sperm, which use a single flagellum to swim; and 
chlamydomonas which uses two fiagella to swim (Bray 
1992). Cilia also occur in very different situations. An 
example is the kinocilium which exists in many hair bun- 
dles of mechanosensitive cells and has the ability to beat 
periodically (Rusch and Thurm 1990). 

The common structural theme of cilia and fiagella is 
the axoneme, a characteristic structure which occurs in 
a large number of very different organisms and cells and 
which appeared early in evolution. The axoneme consists 
of a cylindrical arrangement of 9 doublets of parallel mi- 
crotubules and one pair of microtubules in the center. 




FIG. 1. Schematic representation of the cross-section of the 
axoneme. Nine doublets of microtubules (M) are arranged in 
a cylindrical fashion, two microtubules are located in the cen- 
ter. Dynein motors (D) are attached to microtubule doublets 
and interact with a neighboring doublet. Elastic elements (N) 
such as nexin are indicated as springs. 

How can bending waves be generated by axone- 
mal dyneins which act internally within the axoneme? 
Dynein molecular motors induce relative forces between 
parallel elastic filaments, which are microtubule doublets. 
As a result of these forces, the filaments have the ten- 
dency to slide with respect to each other. If such a slid- 
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ing is permitted globally, the filaments simply separate 
but no bending occurs. Bending results if global sliding 
is suppressed by rigidly connecting the filament pair in 
the region close to one of the two ends. In this situation, 
sliding is still possible locally, however only if the fila- 
ments undergo a bending deformation. This coupling of 
axonemal bending to local microtubule sliding has been 
demonstrated experimentally. In situations where the ax- 
oneme is cut at its basal end, filaments slide and in the 
presence of ATP separate without bending (Warner and 
Mitchell 1981). Small gold-beads specifically attached to 
microtubules in fully functioning flagella can be used to 
directly visualize the local relative sliding during beating 
(Brokaw 1991). 

The theoretical problem of oscillatory axonemal bend- 
ing and wave patterns have been addressed by several 
authors. One can distinguish two principally different 
mechanisms to generate oscillatory forces within the ax- 
oncmc: (i) Deterministic forcing: a chemical oscillator 
could regulate the dynein motors which are activated 
and deactivated periodically. In this case, the regula- 
tory system defines a dynamical force-pattern along the 
axoneme and drives the system in a deterministic way 
(Sugino and Naitoh 1982). (ii) Self-organized beating: 
the axoneme oscillates spontaneously as a result of the 
interplay of force-generating elements and the elastic fil- 
aments (Machin 1963, Brokaw 1975, Brokaw 1985, Lin- 
demann and Kanous 1995). In particular, Brokaw has 
studied thoroughly self-organized patterns of beating 
using numerical simulations of simple models (Brokaw 
1985, Brokaw 1999). Such models are based on the bend- 
ing elasticity of the flagellum and an assumption on the 
coupling of motor activity to the flagellar deformations. 

In the present work, we are mostly interested in self- 
organized beating. Our approach is conceptually differ- 
ent from other works as we focus on an oscillating insta- 
bility of the motor-filament system. In general, sponta- 
neous oscillations occur via a so-called Hopf-bifurcation 
where an initially stable quiescent state becomes unsta- 
ble and starts to oscillate. There is evidence for the 
existence of such a dynamic instability in the axoneme. 
Demembranated flagella show a behavior which depends 
on the ATP-concentration Catp in the solution. For 
small Catp, the flagellum is straight and not moving. If 
Catp is increased, oscillatory motion sets in at a critical 
value of the ATP concentration and persists for larger 
concentrations (Gibbons 1975). This implies an instabil- 
ity of the initial straight state with respect to a wave-like 
mode. Recently, it has been demonstrated using a sim- 
ple model for molecular motors that a large number of 
motors working against an elastic element can generate 
oscillations via a Hopf bifurcation by a generic mecha- 
nism (Julicher and Prost 1995, Jiilicher and Prost 1997). 
This suggests that a system consisting only of molecu- 
lar motors and semiflexible filaments can in general un- 
dergo self-organized oscillations. This idea is supported 



by the facts that flagellar dyneins are capable of gener- 
ating oscillatory motion (Shingyoji et al. 1998) and that 
experiments suggest the existence of dynamic transitions 
in many-motor systems (Riveline et al. 1998, Fujita and 
Ishiwata 1998). Patterns of motion of cilia and flag- 
ella can be complex and embedded in three dimensional 
space. Many examples of propagating bending waves, 
however, are planar and can thus be considered as con- 
fined to two-dimensions (Brokaw 1991). 

In the following sections, we present a systematic study 
of a simplified model for the axoneme which has been in- 
troduced recently (Camalet et al. 1999) and which cap- 
tures the basic physical properties that are relevant for 
its dynamics. In Section II, we present a thorough anal- 
ysis of the dynamics of flexible filaments driven by in- 
ternal forces. This work is inspired by recent studies 
of the dynamics of semiflexible filaments subject to ex- 
ternal forces (Goldstein and Langer 1995, Wiggins and 
Goldstein 1998, Wiggins et al. 1998). We show how the 
dynamic equations can be solved perturbatively and we 
calculate the shapes of bending waves, the velocity of 
swimming and the tension profile along the flagellum. In 
Section III, we briefly review a simple two-state model 
for a large number of coupled molecular motors. This 
model is well suited to represent the dynein molecular 
motors which act within the axoneme. Self-organized 
bending waves and oscillations via a Hopf bifurcation of 
the coupled motor-filament system are studied in section 
IV. We calculate the wave-patterns close to a Hopf bi- 
furcation and determine the frequencies selected by the 
system. Finally, we discuss the relevance of our simple 
model to real axonemal cilia and flagella and propose ex- 
periments which could be performed to test predictions 
that follow from our work. 

II. A SIMPLE MODEL FOR AXONEMAL 
DYNAMICS 

The cylindrical arrangement of microtubule doublets 
within the axoneme can be modeled effectively as an 
elastic rod. Deformations of this rod lead to local slid- 
ing displacements of neighboring microtubules. Here, we 
consider planar deformations. In this case, the geomet- 
ric coupling of bending and sliding can be captured by 
considering two elastic filaments (corresponding to two 
microtubule doublets) arranged in parallel with constant 
separation a along the whole length of the rod. At one 
end, which corresponds to the basal end of an axoneme 
and which we call "head" , the two filaments are rigidly 
attached and not permitted to slide with respect to each 
other. Everywhere else, sliding is possible, see Fig. ||. 
The configurations of the system are described by the 
shape of the filament pair given by the position of the 
neutral line r(s) as a function of the arclength s, where 
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r is a point in two dimensional space. The shapes of the 
two filaments are then given by 



Ti(s) — r(s) — a n(s)/2 
r 2 (s) = r(a) +a n(s)/2 



(1) 



where n with n 2 = 1 is the filament normal. The local 
geometry along the filament pair is characterized by the 
relations 



r = t 
t = C*n 
ri = -Ct 



(2) 
(3) 
(4) 



where t denotes the normalized tangent vector and C — 
t • n is the local curvature. Throughout this paper dots 
denote derivatives with respect to s, i.e. r = d s r = 
dr/ds. 




FIG. 2. Two filaments (solid lines) n and r2 at constant 
separation a are rigidly connected at the bottom end with 
s = 0, where s is the arclength of the neutral line r (dashed). 
Internal forces f(s) are exerted in opposite directions, tan- 
gential to the filaments. The sliding displacement A at the 
tail is indicated. 



A. Bending and sliding of a filament pair 

The two filaments are assumed to be incompressible 
and rigidly attached to each other at one end where s = 0. 
Bending of the filament pair and local sliding displace- 
ments are then coupled by a geometric constraint. The 
sliding displacement at position s along the neutral line 



A(a) 



N) 



(5) 



is defined to be the difference of the total arclcngths along 
the two filaments up to the points ri(s) and V2(s) which 
face each other along the neutral line. From Eqns. ([!]) 
and (0) follows that = (1 ± aC/2)t and thus 



Therefore, 

A = a [ ds'C (7) 
Jo 

is given by the integrated curvature along the filament. 

B. Enthalpy functional 

The bending elasticity of filaments (microtubules) 
characterizes the energetics of the filament pair. In ad- 
dition to filament bending, we also have to take into ac- 
count the large number of passive and active elements 
(e.g. nexin links and dynein molecular motors) which 
give rise to relative forces between neighboring micro- 
tubule doublets. These internal forces can be charac- 
terized by a coarse-grained description defining the force 
per unit length f(s) acting at position s in opposite direc- 
tions on the two microtubules, see Fig. || (internal forces 
must balance). This force density is assumed to arise as 
the sum of active and passive forces generated by a large 
number of proteins. This internal force density corre- 
sponds to a shear stress within the flagellum which tends 
to slide the two filaments with respect to each other. A 
static configuration of a filament pair of length L can 
thus be characterized by the enthalpy functional 



G = 



-C 2 + fA- 



A. 



ds 



(8) 



Here, re denotes the total bending rigidity of the fila- 
ments. The incompressibility of the system is taken into 
account by the Lagrange multiplier function A (a) which is 
used to enforce the constraint r 2 = 1 which ensures that 
s is the arclength (Goldstein and Langer 1995). The in- 
ternal force density / couples to the sliding displacement 
A as described by the contribution J dsfA. The varia- 
tion of this term under small deformations of the filament 
shape represents the work performed by internal stresses. 
Using Eq. M) leads after a partial integration to 



G= / ds 
Jo 



1c 2 -aCF+-r 2 
2 2 



where 



F(s) 



ds'f 



(9) 



(10) 



1*1,2 1 = l±aC/2 



(0) 



is the force density integrated to the tail. From Eq.(|^) 
it follows that if the internal stresses or F are imposed, 
G is minimized for a filament curvature C — Cq where 
Co(s) = aF(s)/n is a local spontaneous curvature. The 
internal forces therefore induce filament bending. In or- 
der to derive the filament dynamics, we determine the 
variation SG with respect to variations Sr. Details of 
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this calculation are given in Appendix A. As a result, we 
find 



Here, 



— = d s [(KC-af)n-Tt] 
or 



t = A + kC — aCF 



(11) 



(12) 



plays the role of the physical tension. This becomes ap- 
parent since from Eq. (Thy) it follows that 



-(«) = t(s) • ^ ds'SG/Sr + ¥ ext {L)^j 



(13) 



where F ext (L) is the external force applied at the end 
which satisfies Eq. (^p|). Therefore, r is the tangent 
component of the integrated forces acting on the filament. 



C. Dynamic equations 

We derive the dynamic equation with the simplifying 
assumption that the hydrodynamics of the surrounding 
fluid can be described by two local friction coefficients £y 
and £j_ for tangential and normal motion, respectively. 
The equations of motion in this case are given by (Wig- 
gins and Goldstein 1998, Wiggins et al. 1998) 



/ 1 1 \ SG 

d t r = - — nn + — tt • — 

U-L f|| J Sr 



(14) 



For the following, it is useful to introduce a coordinate 
system r = (X, Y) and the angle ip between the tangent 
t = (cos ijj, sinip) and the X-axis which satisfies C = tp . 
We find with Eq. Q 



(15) 



d t r = — n ( — k ip +af + tpr) 
+ — t (Kipip — aipf + t) . 



Noting that dtr = ndtip, we obtain an equation of motion 
for ip(s) alone: 



dti) = ~ ( ft 'ip +af + ipi + Tt/j) 
+ —ip(Kipil) - a ftp + t) . 



(16) 



The tension r is determined by the constraint of in- 
compressibility dtr 2 = 2t • d t r = 0. This condition and 
Eq. ([15]) leads to a differential equation for the tension 
profile: 



Eqns. (^6|) and jl7| ) determine the filament dynamics. 
The filament shape follows from 



r(s, t) = r(0, t) + / (cos if), smip)ds' , 



(18) 



where r(0,t) can be obtained from Eq. ( |15| ) evaluated at 
s = 0. 



D. Boundary Conditions 

The filament dynamics depends on the imposed bound- 
ary conditions. The variation SG has contributions at 
the boundaries which can be interpreted as externally 
applied forces F cxt and torques T ext acting at the ends, 
see Appendix A. At the head with s = 0, 



Fext = (kC - a/)n - rt 

J^xt = -kC - a ds'f . 
Jo 

Similarly, at the tail for s = L, 

F cxt = (-kC + af)n + rt 
T cx t = kC 



(19) 



(20) 



If constraints on the positions and/or angles are imposed 
at the ends, forces and torques have to be applied to 
satisfy these constraints. In the following, we will discuss 
different boundary conditions as specified in Table [I]: 

Case A is the situation of a filament with clamped 
head, i.e. both the tangent and the position at the head 
are fixed, the tail is free. Case B is a filament with fixed 
head, i.e. the tangent at the head can vary. The situation 
of a swimming sperm corresponds to case C where the 
friction force of a viscous load attached at the head is 
taken into account, the head is otherwise free. As an 
example of a situation with an external force F ex t(-^) 
applied at the tail we consider Case D. For simplicity, we 
assume a force parallel to the (fixed) tangent at the head. 



boundary 
condition 


head s = 









tail s — L 


A 


d t r = 




d t t 


= 


Fext = 


Text 


= 


B 


d t r = 




Text 


= 


Fext = 


Text 


= 


C 


Fext = - 


(d t r 


Text 


= 


Fext = 


Text 


= 


D 


dtr = 




dtt 


= 


Fext + 


Text 


= 



TABLE I. Different boundary conditions studied. (A) 
clamped head, free tail; (B) fixed head, free tail; (C) swim- 
ming flagellum with viscous load f; (D) clamped head, exter- 
nal force applied at the tail. 



(17) 
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E. Small deformations 

In the absence of internal forces /, the filament relaxes 
passively to a straight rod with -0 = 0, i.e. ip constant. 
Without loss of generality, we choose ip = in this state 
which implies that straight filament is parallel to the X- 
axis. Internal stresses f(s) induce deformations of this 
straight conformation. For small internal stresses, we can 
perform a systematic expansion of the filament dynamics 
in powers of the stress amplitude. 

We introduce a dimensionless parameter e which scales 
the amplitudes of the internal stresses, f(s, t) — efi(s,t), 
where f\ is an arbitrary stress distribution. We can now 
solve the dynamic equations (16) and (^) perturbatively 
by writing 



= e^i + e 2 02 + 0(e 3 ) 
t = t + en + e 2 r 2 + 0(e 3 ) 



(21) 



which allows us to determine the coefficients ip n (s,t) and 
T n (s,t). This procedure is described in Appendix B. We 
find that tp2 and T\ always vanish and To = a is a con- 
stant tension which is equal to the A-component of the 
external force applied at the end. Note, that for bound- 
ary conditions A-C, a = 0. The filament shape is thus 
characterized by the behavior of ^i(s,t) which obeys 



-k ipi +aipi + afi 



(22) 



In order to discuss the filament motion in space, we define 
the average velocity of swimming 



1 r 
lim - / dt'dtY 

S->oo t J 



(23) 



which is independent of s. This velocity is different from 
zero only in the case of boundary condition C. We choose 
the A-axis parallel to v and introduce a coordinate sys- 
tem (x, y) which moves with the filament 



(x,y) = (X-vt,Y) , 



(24) 



where v = ±|v| and the sign depends on whether motion 
is towards the positive or negative A-direction, see Fig. 
||. It is useful to introduce the transverse and longitudinal 
deformations h and it, respectively, which satisfy 



(x,y) = (s + u(s), h(s)) 



(25) 



and which vanish for e = 0. The quantities h, u and v 
can be calculated perturbatively in e, see Appendix B. 
To second order in e transverse motion satisfies 



£±d t h = —ndgh + adgh + ad s f 
Longitudinal displacements satisfy 

u(s) = u(0) - \ I {dshfds 1 



(26) 



(27) 



The boundary conditions for h(s) are given in table 0. 

y aY 




FIG. 3. Resting frame (X, Y) and frame (a;, y) moving with 
the filament at velocity v. The tangent t and the normal n are 
indicated, the angle between the X, a;-axis and t is denoted 
0. 



F. Swimming 

The velocity of swimming v can be calculated pertur- 
batively in e. We consider for simplicity periodic internal 
stresses with frequency u> given byEl 



f(s ) t) = f( S ) i " t + f*(s) ( 



(28) 



which after long time leads to periodic filament motion 
h(s, t) = h{s)e lut + h*(s)e- lut (29) 
where h satisfies according to Eq. ( p6| ) 

Kd^h-adfh + ^iLuh = ad s f ■ (30) 



boundary 
condition 


head s 
F = — 


= 






tail s = L 
d 2 s h = 


A 


h = 




dsh 


= 


hid^h — af 


B 


h = 




nd 2 s h = 


aF 


Kd^h — af 


C 


Cdth = 


af — nd s h 


Kdsh = 


aF 


Kdsh — af 


D 


h = Q 




d s h 


= 


Kdlh 












—ad s h = af 



TABLE II. Boundary conditions for small amplitude mo- 
tion. (A) clamped head, free tail; (B) fixed head, free tail; (C) 
swimming flagellum with viscous load £; (D) clamped head, 
external force applied at the tail. 



1 Note, that more general periodic stresses / = ^ n f n e muJt 
can also be considered. Since the equation of motion ( |2^ ) 
is linear, different modes superimpose linearly and we can 
without loss of generality restrict our discussion to a single 
mode n — 1. 
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FIG. 4. Fundamental modes of filament motion 
hi(s) — He~ 1 ^ close to a Hopf bifurcation are displayed for 
an oscillation frequency of about 50Hz and boundary condi- 
tions A-C. In case C, a viscous load £ ~ 5.10~ 8 Ns/m acts at 
the head, (a) Amplitude H in arbitrary units as a function 
of arclength s. (b) Same plot for the gradient of the phase <j>- 

Writing h = He^ 1 ^ where H and 4> denote the am- 
plitude and phase, respectively, filament motion can be 
expressed as 

h(s,t) = H(s)coa(ut- 0(s)) , (31) 

which represents bending waves for which v p = lj /</> can 
be interpreted as the local wave propagation velocity. Ex- 
amples of bending waves for self-organized beating dis- 
cussed in section IV are displayed in Figs. f| and g. For 
a freely oscillating filament with a viscous load of fric- 
tion coefficient £ attached (case C), we find an average 
propulsion velocity v = Vb/(1 + (/£\\L) where 

"''-(r 1 )^'^ • < 32 » 

is the velocity for ( = 0. Eq. (|3^) which is correct 
to second order in e reveals that motion is only possible 
if the filament friction is anisotropic £||/£± 1 an( i if a 
wave is propagating, 0^0. This is consistent with earlier 
work on swimming (Taylor 1951, Purcell 1977, Stone and 
Samuel 1996). For a rod-like filament with £|| < 
swimming motion is opposite to the direction of wave 
propagation. 




FIG. 5. Filament shapes in the (x, y) plane that correspond 
to the modes displayed in Fig. Snapshots taken at differ- 
ent times to illustrate bending waves are shown for boundary 
conditions A-C. The arrows indicate the direction of wave 
propagation. 

III. INTERNAL FORCES GENERATED BY 
MOLECULAR MOTORS 

The filament dynamics is driven by internal stresses 
f(s,t) generated by a large number of dynein motors and 
passive elements. We describe the properties of a force- 
generating system using a simplified two-state model. 

A. Two state model for molecular motors 

We briefly review the two-state model for a large num- 
ber of motor molecules attached to a filament which slides 
with respect to a second filament (Jiilicher and Prost 
1995, Jiilicher et al. 1997). Each motor is assumed to 
have two different chemical states, a strongly bound state 
1 and a weakly bound state or detached state 2. The 
interaction between a motor and a filament in states 1 
and 2 is characterized by energy landscapes W\(x) and 
Wi (x) , where x denotes the position of a motor along the 
filament. The potentials reflect the symmetry of the fila- 
ments: they are periodic with period I, Wi(x) = Wi(x+l) 
and are spatially asymmetric, Wi(x) ^ Wi{— x). In the 
presence of ATP which is the chemical fuel that drives 
the system, the motors are assumed to undergo transi- 
tions between states. The corresponding transition rates 
are denoted lj% for detachments and loi for attachments. 
We introduce the relative position £ of a motor with re- 
spect to the potential period where x = £ + nl,0<£ i <l 
and n is an integer. The probability to find a motor in 
state i at position £ at time t is denoted Pi(£,, t), P\ + P2 
is normalized within one period. The dynamic equations 
of the system are given by 



G 



d t P\ + vd^Pi = —u\Pi + uj 2 P 2 

d t P 2 + vd £ P 2 = ojiPx - uj 2 P 2 . (33) 

The sliding velocity v = dtx is determined by the force- 
balance 

/ = Xv + p f (Fx dtWi + P 2 d ( W 2 )d( + Kx . (34) 
Jo 

Here, the coefficient A describes the total friction per unit 
length in the system. The number density of motors 
along the filament is denoted by p and / is the force 
per unit length generated by the system. The elastic 
modulus per unit length K occurs in presence of elastic 
elements such as nexins. If motors have an incommensu- 
rate arrangement compared to the filament periodicity, 
Pi + P 2 = l/l and the equations simplify and can be 
expressed by P\ alone 

d t Pi + vd^Pi = -(0J1 + uj 2 )Pi + uj 2 jl 

f = \v + p[ Pid e AW + Kx (35) 
Jo 

where AW = W l -W 2 . 

B. Molecular motors coupled to a filament pair 

The two state model for a large number of motors is 
well suited to represent the internal forces acting within 
the axoneme. We assume that at any position s an inde- 
pendent two-state model described by Eq. (|35| ) is located 
which generates the internal force density f(s). The fil- 
ament sliding displacement is identified with the motor 
displacement, A = x. Note, that here we neglect for sim- 
plicity fluctuations which arise from the chemical activity 
of a finite number of force generators and we assume ho- 
mogeneity of all properties along the axonemal length. 
The energy source of the active system is the chemical 
activity characterized by the transition rates oj\ and uj 2i 
the generated forces induce bending deformations of the 
filaments. 

The dynamic equations ([55]) represent a nonlinear ac- 
tive system which generates time-dependent forces. Since 
we will study periodic motion, we can express forces and 
displacements by a Fourier expansion 

/(f)=Z> eW ( 36 ) 

n 

A(t) = J2 A„e mwt . (37) 

n 

The relation between force amplitudes /„ and amplitudes 
of sliding displacements A„ can in general be expressed 
in powers of the amplitudes A„ (Jiilicher and Prost 1997) 

f n = F^A k+ F^ l A k A l + 0(Al) . (38) 



Here, the summation over common indices is implied. 
The coefficients F^ are calculated for the two-state 
model in Appendix C. However, Eq. (|3|) is more gen- 
eral and characterizes a whole class of active nonlinear 
mechanical systems. 

C. Symmetry considerations 

The interaction of the motors with the two filaments 
is asymmetric. Motors are rigidly connected to one of 
the filaments and slide along the second. In addition, 
the filaments are polar and filament sliding is induced by 
motors towards one particular end of the two filaments. 
As a consequence, the system has a natural tendency to 
create bending deformations with one particular sign of 
the curvature. Exchanging the role of the two filaments 
reverses this sign of bending deformations. This broken 
symmetry does however not reflect the symmetry of the 
axoneme which is symmetric with respect to all micro- 
tubule doublets. Each of the nine microtubule doublet 
within the cylindrical arrangement of the axoneme play 
identical roles and interact with rigidly attached motors 
as well as with motors which slide along their surfaces. If 
all motors generate a constant force, there is consequently 
no resulting bending deformation since all bending mo- 
ments cancel by symmetry. 

(a). 




(b) 

pnnn::: 

FIG. 6. Asymmetric (a) and symmetric (b) mo- 
tors/filament pair. The arrows indicate the polarity of the 
filaments. In case (a), spontaneous bending occurs. In case 
(b) both filaments play identical roles, no spontaneous bend- 
ing occurs. 

In order to introduce this axonemal symmetry in our 
two-dimensional model, we assume that both filaments 
have motors which are rigidly attached and which slide 
on the second filament, see Fig. I The consequence 
of this symmetry is that the filament polarity becomes 
unimportant. Exchanging the role of the two filaments 
is equivalent to replacing the sliding displacement A of 
motors by — A. Therefore, in the context of the two-state 
model this symmetrization is equivalent to the require- 
ment that energy landscapes Wi(x) and transition rates 
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LOi(x) are symmetric functions with respect to x — > —x. 
As a consequence of this symmetry, there is no preferred 
direction of bending. In the following, we adopt this 
choice of a symmetric motor/filament pair. If we use the 
expansion given in Eq. (pq), this symmetry implies that 
all even coefficients F<- 2n >~= 0. 



chemical rates of the ATP hydrolysis cycle. In an exper- 
imental situation, it could be varied e.g. by changing the 
ATP concentration. If the molecular motors are regu- 
lated by some other ion concentration such as e.g. Ca 2+ , 
this concentration could also play the role of the control 
parameter. 



D. Linear response function 



A. Generic aspects 



The perturbative treatment to study small filament de- 
formations introduced in section II can be naturally ex- 
tended to the coupled motor-filament situation using the 
expansion (|3^). Up to second order in e, only the lin- 
ear response coefficient of the active system plays a role 
which for the two-state model is given by FjQ = x$nk 
with 



X = K + \iuj- p [ d£ deAW diR , 
Jo ' a + iuj 



(39) 



Here, we have introduced R = Lo 2 /al and a — u)% + lu 2 . 
Higher order terms F^ 2n+1>> have to be taken into account 
if the third or higher order in e is considered. The linear 
response function x as wen as nonlinear coefficients can 
be calculated most easily for a simple choice of symmetric 
potential and transition rates 



AW(£) = [/cos(2tt£/Z) 

=j8-0cos(2tt£/O 
=a-j0 + /3cos(27rf/Z) 



(40) 



where a and /3 are ^-independent rate constants. For this 
convenient choice 



x(Sl, u>) = K + i\lo — pkQ 



iuj/a + [oj/aY 
1 + (cu/a) 2 



(41) 



where k = U/l 2 is the cross-bridge elasticity of the mo- 
tors. We have introduced the dimensionless parameter 
SI = 2n 2 (3/a with < fl < ir 2 which plays the role of 
a control parameter of the motor-filament system, a is a 
characteristic ATP cycling rate. 



IV. SELF-ORGANIZED BEATING VIA A HOPF 
BIFURCATION 

A Hopf-bifurcation is an oscillating instability of an ini- 
tial non-oscillating state which occurs for a critical value 
Sl c of a control parameter. For SI < Q c , the system is 
passive and not moving, while for SI > Sl c it exhibits 
spontaneous oscillations. As we demonstrate below, self- 
organized oscillations of the driven filament pair are a 
natural consequence of its physical properties. The con- 
trol parameter SI is in our two-state model a ratio of 



For oscillations in the vicinity of a Hopf-bifurcation, 
filament deformations are small and we can thus use Eq. 
( pf ) to describe the filament dynamics. The force ampli- 
tude / can be obtained from the expansion 



;=X(S1, W )A + 0(A 3 ) 



(42) 



in the amplitude of sliding displacements. In addition, 
this amplitude is related to the filament shape: 



A(s) = a(d s h(s) - d 3 h(0)) + 0{h 3 



(43) 



With Eq. (poj) , we find that spontaneously oscillating 
modes h(s) at frequency w are to linear order near the 
bifurcation solutions to 



ndjh - ad 2 h + ^iwh = a 2 X d 2 h 



(44) 



Note that this equation is general and its structure does 
not depend on the specific model chosen for the force 
generating elements. Boundary conditions corresponding 
to cases A-D follow from Table [n] by replacing h —> h, 
dth — > ibjh and / — > / = ax{d s h — d s h(0)). 

As discussed in Appendix D, Eq. (|44|) has the struc- 
ture of an eigenvalue problem, with \ playing the role 
of a complex eigenvalue. For every choice of parameters, 
there exists an infinite set of nontrivial solutions h n to 
Eq. (44) if x is equal to the corresponding eigenvalue 
Xn, n — 1,2, ... We order these values according to their 
amplitude: \x n \ < |Xn+i|- The functional dependence 
of the Eigenvalues Xn on the model parameters is com- 
pletely determined by dimensionless functions Xn (^, u>) 
according to 



Xn 



L 2 a 2 



Xn{(7,Lu) 



(45) 



where a = <jL 2 /k and u> = u>£,±L 4 /k are a dimensionless 
tension and frequency, respectively. 

The existence of discrete eigenvalues reveals that only 
particular values of the linear response function of the 
active material are consistent with being in the vicinity 
of a Hopf-bifurcation. If the actual value of x(S1,cj) of 
the system differs from one of the eigenvalues Xn, the 
system is either not oscillating, or it oscillates with large 
amplitude for which higher order terms in e cannot be 
neglected. An important consequence of this observa- 
tion is that the modes h n of filament beating close to a 
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Hopf-bifurcation can be calculated without knowledge of 
the properties of the active elements. Not even knowl- 
edge of the linear response coefficient is required since it 
follows as an eigenvalue at the bifurcation. Therefore, 
filament motion close to a Hopf bifurcation has generic 
or universal features which do not depend on details of 
the problem such as the structural complexity within the 
axoneme. The motion is given by one of the modes h n 
which only depend on a and u). 







Imx 
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FIG. 7. (a) Eigenvalues \ n for boundary conditions B. 
Each eigenvalue is represented by a line which traces the loca- 
tions of x» in the complex plane for varying reduced frequency 
Q. The lines begin for uo = on the real axis, the value of ui 
at the other ends are indicated. 

Fig. |I] displays the amplitude H and the gradient of 
the phase of the fundamental mode hi = He~ lc ^. Note, 
that we use arbitrary units for H since the amplitude of 
h(s) is not determined by Eq. (f44|). Corresponding fil- 
ament shapes represented according to Eq. ( poj ) in the 
(x, y) plane, are displayed in Fig. |j| for boundary con- 
ditions A-C as snapshots taken at different times. The 
oscillation amplitudes are chosen in such a way that the 
maximal value of ip ~ n/2. The direction of wave propa- 
gation depends on the boundary conditions. For clamped 
head (case A) , waves propagate towards the head whereas 
for fixed or free head (cases B and C), they propagate 
towards the tail. The positions of the lowest eigenvalues 
Xn hi the complex plane for boundary conditions B are 
displayed in Fig. [7] for a = and varying uj. 



B. Selection of eigenmodes and frequency 

In order to determine which of the modes h n is selected 
and at what frequency oj the system oscillates at the bi- 
furcation point, explicit knowledge of the linear response 
function x(^,^) is necessary. As a simple example, we 
use x as given by Eq. ( |4l"| ) for a two-state model. For 

= 0, x = K + iXoj which is a passive viscoelastic re- 
sponse. In this case, no spontaneous motion is possible 



which can be seen by the fact that all eigenvalues %n 
have negative real and imaginary parts which cannot be 
matched by the linear response of a passive system. If Q 
is increased, we are interested in a critical point fi = Q c 
for which the straight filament configuration becomes un- 
stable. This happens as soon as the linear response func- 
tion x matches for a particular frequency u> — lo c one of 
the complex eigenvalues, 



a 2 L 2 



Xn(cr,uj c ) = x(ttc,u c ) 



(46) 



where uj c — £j_L 4 cj c /k. Since x is complex, both the real 
and the imaginary part of Eq. (M) represent indepen- 



dent conditions. Therefore, Eq. ( |46| ) determines both 
the critical point fl c as well as the selected frequency io c . 
The selected mode n is the first mode (beginning from 
Q = 0) to satisfy this condition. Since |x(fi,w)| typi- 
cally increases with increasing f2, the instability occurs 
almost exclusively for n = 1 since |xi| has been defined 
to have the smallest value. For f2 > f2 c , but close to the 
transition, the system oscillates spontaneously with fre- 
quency uj c . The shapes of filament beating are character- 
ized by the corresponding mode h n (s) whose amplitude 
is determined by nonlinear terms. For the case of a con- 
tinuous transition, the deformation amplitude increases 
as \h n \ ~ {n-VLc) 1 ' 2 . 



C. Axonemal vibrations for different lengths 

We choose the parameters of our model to correspond 
to the axonemal structure. We estimate the bending 
rigidity n ~ 4 • 10~ 22 Nm 2 which is the rigidity of about 
20 microtubules. Furthermore, we choose a microtubule 
separation a ~ 20nm, motor density p ~ 5 ■ 10 8 m _1 , 
friction per unit length A ~ lkg/ms, rate constant 
a ~ 10 3 s _1 , cross-bridge elasticity k ~ 10 _3 N/m and 
a perpendicular friction £j_ ~ 10~ 3 kg/ms which is of 
the order of the viscosity of water. A rough estimate 
of the elastic modulus per unit length associated with 
filament sliding can be obtained by comparing the num- 
ber of dynein heads to the number of nexin links within 
the axoneme. This suggests that K is relatively small, 
K<kp/10. 

The selected frequency lu c at the bifurcation point for 
case B is shown in Fig. |^ as a function of the axoneme 
length for K = 0. For small lengths, the oscillation fre- 
quency is large and increases as L decreases. In this 
high-frequency regime, which occurs for L^10/im, the 
system vibrates in a mode with no apparent wave prop- 
agation. For L^10/mi the frequency is only weakly In- 
dependent and the system propagates bending waves of 
a wave-length shorter than the filament length. 

The limit of small lengths corresponds to small Co and 
can be studied analytically. For a = 0, the eigenvalue 



Xi is given to linear order in w by xi 



-7T /4 — ijuj. 
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The coefficient 7 depends on boundary conditions but 
not on any model parameters. For a clamped head (A), 
7 ~ 0.184, for a fixed head (B), 7 ~ 0.008, see Appendix 
D. The criterion for a Hopf bifurcation for small L is 



X(^c,^c) 

Together with Eq. 



Aa 2 L 2 



(47) 



J), we find the critical frequency 



7T 

21? 



na 
7£i 



1/2 



1 + (L/L K f 



1 



1/2 



(48) 



(W£) 2 . 

where L K = (n 2 n/AKa 2 ) 1 / 2 and L x = (Xa 2 /jC±) 1/2 are 
two characteristic lengths. For K and A small, Lk 
L ^> L\ and the critical frequency behaves as lo c ~ 1/L 2 . 




L(um) 



FIG. 8. Oscillation frequency u c /2tt at the bifurcation 
point as a function of the axoneme length L for boundary 
conditions B and model parameters as given in the text. The 
insets show characteristic patterns of motion for small and 
large lengths. 

The critical value of the control-parameter for small L 
is given by 

n c ~ ^ K + Xa 1 a ^ L " . (49) 
Aa 2 pkL 2 pk pka 2 

For axoneme lengths below a characteristic value, L < 
L m j n , the condition Q, c < ir 2 is violated. Therefore, os- 
cillations exist only for L > L mul . If (A" + Xa)/pk <C 1, 
L m i n ~ yj K/4a 2 pk ~ 1/im for our choice of the param- 
eters. The maximal frequency of oscillations occurs for 
L = -Lmin- Within the range of lengths between lpm and 
50pm, we find frequencies lu c /2tt which vary between tens 
of Hz and about 20kHz. The parameters chosen in Fig. 
H lead to frequencies above 150-ffz. Lower frequencies 
are obtained for larger values of the friction A or smaller 
chemical rate a. 



D. Effect of external forces applied at the tail 

In the presence of an external force applied at the end 
a / and the filament pair is under tension. The eigen- 



values Xn depend on a, therefore the tension affects shape 
and frequency of an unstable mode. Inspection of Eq. 
(H) suggests that the tension a can be interpreted as an 
additive contribution to the linear response coefficient 
X — » X + a /a 2 . Taking into account the boundary condi- 
tions complicates the situation, the tr-dependence of the 
eigenvalues is in general nontrivial. However, for a fila- 
ment with clamped head and a force applied at the tail 
parallel to the tangent at the head (case D) the effect of 
tension can be easily studied. In this particular case, the 
tension leads to the same contribution to x both in Eq. 
( p4| ) and in the boundary conditions and the eigenvalues 
thus depend linearly on a: 



X„(cr,w) = Xn(0,U)) 



(50) 



Since the tension corresponds to a change of the real part 
of x, its presence has the same effect on filament motion 
as an increase of the elastic modulus K per unit length in 
Eq. ( fi~i|) by a /a 2 . Therefore, we can include the effects 
of tension on the critical frequency in Eq. ( fi"8|) by replac- 
ing K — » K + a /a 2 which reveals that in those regimes 
where the frequency depends on A, it will increase for 
increasing tension. The critical value of the control pa- 
rameter Q c is a function of the tension, with Eq. (fiT 
we find 



n c (a) = n c (o) 



1 



pka 2 



(51) 



For our choice of parameters, pka 2 ~ 200pN which in- 
dicates that forces of the order of 10 2 pN should have a 
significant effect on the bifurcation. Furthermore, Eq. 
( pl| ) indicates that the tension a can play the role of a 
second control parameter for the bifurcation. Consider 
a tensionless filament oscillating close to the bifurcation 
for f2 > f2 c (0). If a tension is applied, the critical value 
f2 c (<7) increases until the system reaches for a critical 
tension a c a bifurcation point with Q c (a c ) — and the 
system stops oscillating. 



V. DISCUSSION 

In the previous sections, we have shown that many 
aspects of axonemal eating can be described by a simple 
model based on the idea of local sliding of microtubule 
doublets driven by molecular motors inside the axoneme. 
This model represents a class of physical systems termed 
internally driven filaments (Camalet et al. 1999) which 
have characteristic properties that are closely linked to 
the geometric constraint that couples global bending and 
local filament sliding. 

Our work shows that axonemal beating can be studied 
both numerically and analytically using a coarse-grained 
description which ignores many details of the proteins in- 
volved and is based on effective material properties such 
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as bending rigidities, internal friction coefficient and elas- 
tic modulus per unit length as well as the frequency de- 
pendent linear response function of the active elements. 
Our main results are: (i) The filament pair introduces a 
geometric coupling of the active elements at different po- 
sition along the filaments. This coupling is suited for the 
generation of periodic beating motion by a general self- 
organization mechanism of the system, (ii) The genera- 
tion of oscillations and propagating bending waves occurs 
via a Hopf bifurcation for a critical value of a control- 
parameter such as the ATP concentration. The patterns 
of motion generated close to this bifurcation can be cal- 
culated without any knowledge of the internal force gen- 
crating mechanism if the oscillation frequency is known, 
(iii) The frequency depends on chemical rates of the mo- 
tors, coefficients of internal elasticities and frictions, the 
solvent viscosity and the microtubule rigidity. Our cal- 
culations using a simple two-state model for molecular 
motors suggest significant variations of the oscillation fre- 
quency if the axoncmal length is varied. Short axonemes 
are predicted to be able to oscillate at frequencies of sev- 
eral kHz. This fact is particularly interesting in the case 
of kinocilia which are located in the hair bundles of many 
mechanosensitive cells that are involved in the detection 
of sounds. If cilia are able to vibrate at high frequencies, 
this suggests that the kinocilium could play as active role 
in sound detection (Camalet et al. 2000) . (iv) Forces ap- 
plied to the axoneme can be an important tool to study 
the mechanism of force generation. An external force 
can play the role of a second control parameter for the 
bifurcation which has a stabilizing effect for increasing 
forces. 

These results have been obtained using several simpli- 
fying assumptions. Our model represents the solvent hy- 
drodynamics by an anisotropic local friction acting on the 
rod- like filament. This approach ignores hydrodynamic 
interactions between different parts of the filament which 
lead to logarithmic corrections. These hydrodynamic ef- 
fects do not change the basic physics but can lead to cor- 
rections to the numerical results. It is straightforward to 
generalize our model to incorporate the effects of hydro- 
dynamic interaction, but this does not change our results 
qualitatively. A more serious simplification of our model 
is the restriction to a two-dimensional system and to fila- 
ment configurations which are planar. This choice is mo- 
tivated by the fact that observed bending waves of flagella 
are planar in many cases (Brokaw 1991). If filament con- 
figurations arc planar then our model applies and is com- 
plete. However, this model cannot explain why motion 
is confined to a plane and it misses all non-planar modes 
of beating. In order to address such questions, a three 
dimensional generalizations have to be used. Such gen- 
eralizations introduce additional aspects. In particular, 
torsional deformations become relevant. The local sliding 
displacements depend in the three dimensional case both 
on bending and torsional deformations and the full tor- 



sional dynamics has to be accounted for in the dynamic 
equations. Finally, we have restricted most calculations 
to the limit of small deformations which corresponds to 
filament shapes that are almost straight. This regime has 
several important features. The filament dynamics can 
be studied by an analytic approach by a systematic ex- 
pansion in the deformation amplitude. This allows us to 
characterize linear and nonlinear terms both for the fila- 
ment dynamics as well as for the properties of the active 
elements. Patterns of motion close to a Hopf bifurcation 
are fully characterized by the linear terms of this expan- 
sion. These linear terms are given by the structure of the 
problem and their form does not depend on molecular 
details of the system. For example, most details of the 
operation of the molecular motors are unimportant for 
the shape of the filament oscillations at the bifurcation, 
only the linear response function of the active material 
plays a role. 

If filament beating with larger amplitude is of interest, 
nonlincarities become relevant. Nonlincarities arise due 
to nonlinear geometric terms in the bending energy or via 
non-Hookian corrections of the elasticity. Furthermore, 
nonlincarities in the force-generation process of molecular 
motors exist. All these nonlinearities determine the large 
amplitude motion and could give rise to new types of be- 
havior such as additional dynamic instabilities. However, 
in contrast to the linear terms, the form of the dominant 
nonlincarities does depend on structural details of the ax- 
oneme. Therefore, an analysis of large amplitude motion 
is difficult and knowledge of the nature of the dominating 
nonlincarities is required. However, if no new instabili- 
ties after the initial Hopf-bifurcation occur, the princi- 
pal effect of nonlinear terms is to fix the amplitude of 
propagating waves. We therefore expect that propagat- 
ing waves with larger amplitude are in many cases well 
approximated by our calculation to linear order. 

Our work shows that propagating bending waves and 
oscillatory motion of internally driven filaments can oc- 
cur naturally by a simple physical mechanism. Complex 
biochemical networks to control the system are not re- 
quired for wave propagation to occur. This suggests that 
the basic axonemal structure intrinsically has the ability 
to oscillate. Biochemical regulation systems are thus ex- 
pected to control this activity on a higher level but are 
not responsible for oscillations. Our work shows that ex- 
perimental studies of axonemal beating close to a Hopf 
bifurcation would be very valuable. Such experiments 
could be e.g. performed by using demembranated flag- 
ella and ATP concentrations not far from the level where 
beating sets in. The observed behavior close to the bifur- 
cation would give insight in the self-organization at work. 
Furthermore, externally applied forces and manipulation 
of boundary conditions could be sensitive tools to test 
some of the predictions of our work. 
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APPENDIX A: FUNCTIONAL DERIVATIVE OF 
THE ENTHALPY 

The variation SG of the enthalpy G given by Eq. (||) 
under variations Sr of a shape r(s) is 



5G= [ [{kC - aF) 6C + Ar • Sr ] ds 
Jo 



(Al) 



Note, that f 2 = 1 but (r+<5r) 2 1 in general. Under such 
a variation, SC — <5(n-r) = n-Sr since r-<5n = Cn dn = 0. 
Therefore, 



SG = [ [(kC -aF)n-Sr + At-Sr] ds 
Jo 

Two subsequent partial integrations lead to 
SG = [(kC -aF)n- Sr}^ 
+ (-(kC - af)n + {kC 2 - aCF + A)t) • Si 



(A2) 



(A3) 



+ / d s ((KC-af)n- (kC 2 -aCF + A)t) • Sr ds 



The functional derivative SG/Sr given by Eq. ( |ll| ) can be 
read off from the integrand, the boundary terms provide 
expressions for forces and torques at the ends given by 
Eqs. ph and Q. 



APPENDIX B: SMALL DEFORMATIONS 

Inserting the expansions ( pl| ) in the differential Eq. 
( |l7| ) for the tension profile leads at each order in e to a 
separate equation. Up to second order we find 



TO =0 

n = o 



(Bl) 



f-2 = d s (ipi(afi - Ktpi)) + —ipi(-K -0i +afi + T ipi) ■ 

SJL 

From the boundary conditions, it follows that tq = a is 
constant and n = 0. Repeating the same procedure for 
the dynamic Eq. (|16|) using n = and tq = <r, we obtain 



£,±d t *l>2 = -K ^2 



(B2) 



The Equation for tp2 is independent of f\ and only de- 
scribes transient behavior which depends on initial con- 
ditions. After long times and for limit cycle motion, 



Tp2 = 0. The transverse and longitudinal displacements 
h, u and the velocity v are obtained perturbatively as 



h = ehi + ( 2 h 2 + 0(e 3 ) 
u = eui + ( 2 u 2 + 0(e 3 ) 
v = evi + e 2 v 2 + 0(e 3 ) 



Using Eqns. (21) and ( |lq ) we find ui(s) 
/i2(0) and 



h 1 (s) = h 1 (0) 



ipi{s')ds' 



1 



u 2 («) = ua(0) - | ^ Ms'fds' 



(B3) 

ui(0), /i 2 (s) = 

(B4) 
(B5) 



The dynamics of h(s) and u(s) to second order in e is 
therefore determined by ipi(s,t) and the motion <9 t r for 
s = 0. For the latter, we find from Eq. ( |l5| ) «i = 0, 
StUi(O) = 0, 9 t /i 2 (0) = and 

e±ftfci(0) = -« )ii (0) + T hi(0) + o/i(0) (B6) 
f||(«a +flt«a(0)) - -e||0i(O)9f/ii(O) + ^1(0)^(0) 

- 0^(0)^(0) + Ta(0) . (B7) 

For the lateral motion, we obtain with Eq. flB2| ) and 



£±(dthi(s) - d t hi(0)) = -k hi +ahi + afi 



Using Eq. (B6), we obtain 

£±dthi = — k h\ +ahi + afi 



(B8) 



(B9) 



which is the result given in Eq. ( |2q) . In order to calculate 
t 2 (s) and V2 we integrate Eq. (Bl). Together with Eq. 



T2(s)=T 2 (0)+ "0l(a/l - K?/>l) 

+ / Ad t hids' . 



(B10) 



After elimination of 72(0) via Eq. (B7) we find the 
tension profile 

t 2 (s) = T 2 (0) + s^(v 2 + dtu 2 (0)) 



ds' 



ipi(afi - n-ipi) + £\\ipid t hi 



£|l / ds' / ds"i/jidttpi 



(Bll) 



The tension at the head 12 (0), as well as v 2 + dtu 2 (0) are 
fixed using the boundary conditions. In cases A,B and 
D where the head is fixed V2 + <9t 1*2(0) = 0. In this case, 
72(0) follows from Eq. (Bll) at s = L. For boundary 
condition C, the velocity of the head is determined from 
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the condition T2(L) = together with the boundary con- 
ditions r 2 (0) - -Vi(0)M 1 (0)- a / 1 (0)) + C(v 2 +5 t u 2 (0)) 
and k hi (L) = afi(L) which leads to 

v 2 + d t u 2 {0) = t X .~Jl [ htdthxds 



(D2) 



1=1 



In order to determine the amplitudes A the four bound- 
ary conditions have to be used. Since Eq. ( fli]) is linear 
and homogeneous, the equations for the coefficients A, 
have the form 



ds / ds'-d t {h\) 



(B12) 



Averaging this equation for periodic motion h\{s,t) 
H(s) cos(ujt - <p(s)) leads to Eq. Q. 



AiMij = 



(D3) 



APPENDIX C: LINEAR AND NONLINEAR 
RESPONSE FUNCTIONS OF THE TWO-STATE 
MODEL 

We introduce the Fourier modes of the distribution 
function P x (£, t) = J2 k p i (6 fc)e lfcw * . The dynamic equa- 
tion (HH) of the two-state model can then be expressed 
as 

lS k 



where My is a 4 x 4 matrix which depends on x, & ancl 
Nontrivial solutions exist only for those values x = Xn 
for which det My = 0. The c orre sponding eigenmode 



is the solution for Aj of Eq. (D3). These eigenvalues 
and eigenfunctions can be determined by numerically ob- 
taining solutions for det My = 0. In the limit of small 
ui, analytic expressions can be obtained by inserting the 



ansatz x = X 



(o). 



■ X & ' 



0{lo ) for the eigenvalue which 
-,i/2 



leads to an expansion of det My in powers of u> ' . This 

procedure leads for a = to Xn = — [i n — l) 71 " + 7r /2] 2 
independent of boundary conditions. The linear coeffi- 

(i) _ 



Pi(£, k) = 5 k . R(n - iuV k,l t , AidtPifc, m) (CI) cient Xi = ~*7 depends on boundary conditions. For 



Inserting the Ansatz 



boundary conditions A 

12tt 3 - 32tt 2 + 1728tt - 4608 



Pr(£, fc) - ay fci0 + Pij ) A l + PiZ^m + 0(A 3 ) (C2) 



7T 5 + 1447T 3 



0.184 



(D4) 



while for case B, 



into Eq. (CI), one obtains a recursion relation 



lSk,k n +l R p(n-l) 



192 /tt 
7 = — 



32 12 



tt j + 2tt^ - 2tt ~ 0.008 . (D5) 



(C3) 



(n) 

the knowledge of the functions Pf, fel fe allows one to 
calculate the coefficients of the expansion (|3^). In par- 
ticular, with P( ' = i?, one obtains the result (39) and 



(n) 



fe.fei ,..,k n 



Jo 



(C4) 



APPENDIX D: EIGENVALUES AND 
EIGENMODES NEAR A HOPF BIFURCATION 

Nontrivial solutions to Eq. (|4|) can be found using the 
Ansatz h = Ae qs l L with q being solution to 



(ff + X)9 2 + »w = 



(Dl) 



Here, <7 = <jL 2 /k, uj = to^±L /re, and % = x a2 L 2 /K are 
a dimensionless tension, frequency and linear response 
coefficient, respectively. Eq. (Dl) has four complex so- 
lutions qi, i = 1, .., 4. Therefore, 
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